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The Ginzburg-Landau (GL) theory is recast using a Hamiltonian involving the complete kinetic 
energy density which requires that the surface energy must contain a term V| , )/'| 2 to support super- 
conducting (SC) states. The GL equations contain two temperature, t, dependent parameters a(t) 
and f3(t), which are respectively the coefficients of the SC pair density oc \ip\ 2 , and the pair inter- 
action term oc \i/j\ 4 in the free energy density. The sign of these parameters, which define distinct 
solution classes, and the ratio s(t) = y^a|/|/3| are governed by the characteristics of the surface 
energy density. In addition to the conventional bulk superconducting states with (a < 0, (5 > 0), 
anomalous superconducting states exist for all other sign combinations, including cases with (3 < 
which may exist only when surface pair interactions are significant. All possible solutions of our 
generalized nonlinear, one dimensional GL equations are found analytically and applied to a thin 
superconducting slab which manifests the possibility of states exhibiting enhanced, diminished, and 
pre-wetting superconductivity. Critical currents are determined as functions of s(t) and surface 
parameters. The results are applied to critical current experiments on SNS systems. 



I. INTRODUCTION 

Macroscopic effects of superconductivity are investigated using the Ginzburg-Landau (GL) phenomenological model 1 
for the free energy of the superconducting state characterized by a complex order parameter ip and two parameters 
a and (3. These parameter are, respectively, coefficients of the linear superconducting (SC) pair density term oc \ij)\ 2 
and the quadratic pair interaction term oc l^l 4 - A variational process leads to coupled GL differential equations that 
characterize a minimum free energy state that depends on the surface energy density through additional parameters. 
In almost all cases the source of an applied magnetic field or a current is externally controllable and the free energy 
is of the Gibbs form. Quantized magnetic levitation 2 is an exceptional, fundamental example of a nano-magneto- 
mechanical device, whose analysis requires the use of the Hclmholtz free energy. This article deals with the Gibbs free 
energy and the anomalous superconductivity induced by externally controlled surface energy. Since our Hamiltonian 
functional differs from the standard GL functional, and since unexpected phenomena arising from the surface energy, 
the article is written in a self-contained format with all relations derived from our energy functional. 

In the absence of surface effects, the pair interaction parameter (3 must be positive in order to have an energy 
minimum and a concomitant SC state. The surface boundary condition of the superconducting order parameter 
modifies the free energy density near the surface on the scale of the coherence length. For sample dimensions of Imm 
and larger, boundary conditions do not contribute significantly to the total free energy. However, for a SC sample, 
assumed to be uniform and homogeneous, with relevant dimensions of 10/im and smaller, the total free energy can be 
significantly altered by boundary effects. 

For a weakly superconducting surface, considered previously 3-5 , the surface energy term added to the Ginzburg- 
Landau (GL) free energy, has the form oc J s (\ip\ 2 /b2)ds, where the integral is extended over the sample surface, and 62 
is a characteristic length. When 62 < 0, the overall free energy is lowered and an increase of the transition temperature 
above the bulk value T c is possible 6 . If 62 > 0, the free energy is increased and the nucleation temperature is lowered 
. - - 1 relative to T c . In the weak SC surface limit, the sample is superconducting only if f3 is positive. 

If the surface is strongly superconducting, the surface energy has an additional SC pair interaction term 
oc j s (\ip\ A /b^ds, with 64 a characteristic length. For 64 > 0, we show that the entire sample may be in an anomalous 
superconducting state with (3 negative. We also show that the minimum free energy explicitly depends on 64; whereas 
the dependence on 62 is only implicit through Thus, the parameter 64 essentially determines whether a partic- 
ular solution of the differential GL equations represents a physical SC state, or whether it is merely a mathematical 
solution. 

To account for the spatial variation of ip, the conventional GL energy functional includes the term oc Ip^l 2 , where 
p is the momentum operator. In this study we introduce a Hamiltonian with the complete kinetic energy density 
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K oc ip*p 2 tp = \ptp\ 2 — ieW ■ (ip*pip), where the parameter e = 1, or is used to keep track of the difference between 
the K model and the standard GL model. In view of the divergence operation, the additional term is a surface 
term. The imaginary part of the divergence term is identically zero, but the real part changes the surface boundary 
conditions, and hence the character of the order parameter throughout the sample volume. It is shown that the real 
part of the divergence term in the K model destroys superconductivity unless the surface energy is augmented by a 
term c|r/'|V|i/'|) where c > e is a constant. Furthermore, a second order phase transition is possible only if c = e, and 
the transition temperature is higher for c = e than for any value of c > e. For c < e SC states do not exist for small 
values of the order parameter; however, first order phase transitions may exist. 

In the GL theory a and (3 are arbitrary parameters, assumed to be functions of the reduced temperature t = T/T c . 
However, the explicit temperature dependence of a and (3 is not given and thus has to be conjectured from experiments. 
Since we do not want to assume a priori an explicit temperature dependence, we first investigate the solutions of the 
GL equations consistent with the boundary conditions, independent of specific temperature models for a and (3. For 
one dimensional systems with a uniform current density j, all possible physical solutions of the GL equations, with 
sgn(a) = ±1 and sgn(/3) = ±1 are found, and categorized. These solutions are then applied to a plane slab SC sample 
and the influence of the surface energy parameters 62 and 64 on the sample order parameter is studied, and unusual 
SC states, including "pre- wetting" surface states, are analyzed. Critical currents are shown as functions of s(t). A 
temperature model is introduced for the fundamental ratio s{t) = y/\a{t)\/\j3{t)\ to illustrate possible enhancement or 
reduction of the critical temperature. Finally, the critical current of a superconductor-normal-superconductor (SNS) 
layered system is calculated, and the results with a > 0,(3 > are shown to be in excellent agreement with critical 
current experiments. 7 



II. GENERAL EQUATIONS 



As originally formulated by Ginzburg and Landau 1 , the free energy of a superconductor in a magnetic field H is 
a function of a coordinate dependent complex order parameter tp(r) and the vector potential A(r). We begin our 
development with a general energy functional G that encompasses a broad spectrum of models, including those leading 
to the stationary states of the Schrodinger equation and the GL equations for SC. It is 
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where U is a potential energy density that depends only on functions of r and |i/>(r)|, K(r) is the kinetic energy 
density, and A' is a linear function of A. Specific forms of U and A' will be detailed later. The second integral is 
over the surface S enclosing the volume V. The unit vector n is normal to the surface at each point on S. The surface 
integrand function A, which has units of (length) -4 , is assumed to contain all parameters that characterize not only 
the intrinsic aspects of the surface, but also external effects that link to the sample volume via the surface. It is shown 
below that the cV|V>| 2 term is necessary for a SC state to exist, and the value of c determines whether a SC state 
exists and whether a phase transition to the SC state is first or second order. The boundary conditions for ip depend 
on c and A, and a specific form for A is introduced. 

The kinetic energy density if of a particle of effective mass to* and charge magnitude e* subject to a vector potential 
A is given by 
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K(r) = ip*p 2 ip = |pVT - *eV • {ip*pip), 



(2.2) 



where 



P = -«V + —A, 

n 

The first term in Eq. (2.2) is the form used in the standard GL functional. To distinquish the effects arising from the 
complex divergence term the parameter e = 1 or is introduced. When e = one recovers the standard GL model. 
As shown below, only the real part of the divergence term in K contributes to G. 
Introducing the modulus-phase form 



^(r) = Mr)|exp[*0(r)], 



(2.3) 



the density K assumes the form 
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-^K(r) = [V|V| • VM + |^| 2 Q • Q] - 26V • [|^| 2 Q - -V|V| 2 ], (2.4) 
where the vector Q is the gauge invariant phase gradient 

with </> c = /i/e* the flux quantum. In the superconducting case, e* = 2|e| > is the SC pair charge. The magnetic 
field H is given by 

HoH = V x A = x (Q - V0). (2.6) 

As will be shown, the term V • [|V1 2 Q] = in Eq. (2.4), since it is proportional to the divergence of the real current 
density. The term V • V|-0| 2 may appear to be the divergence of a diffusion current, but it is not. It contributes to 
the surface energy. Using Eq. (2.4), the functional G is 
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/ ds • [nA(|^|,6>, A) + )-{a- 1)V|V>| 2 ], with a=l + c-e. 
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The conventional GL free energy difference AG = G sc — G n between the superconducting and normal states is given 
by G in Eq. (2.7) with c = e = and A' = A — A , where A a is the vector potential of the applied field H a . Note 
that for the standard GL model e = in K oc tp*p 2 tp — \ptp\ 2 — ieW ■ (ip*pip). The vector potential A that appears 
in Q reduces to A a for single particle models, since a single charged particle does not interact with its own field. 

The vectors V0 and A in Q arc distinct mathematical entities. When a magnetic field is present V x A ^ 0, 
whereas V x V(9 = 0, except at vortex centers, where V x V6> is a sum of delta functions 8 ' 9 . Thus 9 and A must be 
considered as independent variables, and the minimum set of independent variables contains three elements, which 
we chose as (|^>|,0, A). As outlined the Appendix, setting the variation SG = 0, leads to the differential equations in 
the volume and corresponding surface boundary conditions. The variation with respect to \tp\,9 and A, respectively 
yield 

*w = i?m +m - Q ' (2 - 8) 

v • (|VI 2 Q) = 0, (2.9) 
-Vx (Vx A') = -6 S ^|^| 2 Q. (2.10) 

fl TYl 

The parameter = for single particle models and e s = 1 for models that include particle interaction terms, e.g. 
the GL model. 

The corresponding surface boundary conditions are 

[(« + l)VM + ( fl -l)M^]-n = -^ (2.11) 

2|^| 2 Q-n=— , (2.12) 

and 

1 h 2 

— n x (V x A') = ^V A A. (2.13) 

As outlined in the last paragraph of the Appendix, the minimum energy functional consistent with Eq. (2.8)- (2.13) 
assumes the form 
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Using the boundary condition (2.11) with a = 1, the expression for G m i„ can be written in the form 

G m in = J^r {g(U) + ^[V x Af } + ^ / & ■ n 5 (A), 



(2.14) 



(2.15) 



where 



Equation (2.15) shows the functional symmetry of the volume potential energy density U and the surface potential 
energy per unit area (h 2 /2m*)A which determine G m i n - It is evident from this form of G m i n that if U and A arc 
expanded in powers of the particle density the coefficient of leading term \ip\ 2 does not contribute directly to G m i n - 
To proceed further one must introduce models for the functions U and A. 
For the potential energy density U we assume the form 



t/(r,|^(r)|) =a(r)|</>(r)| s 



l/3(r) 



(2-16) 



where n* is a constant reference particle density. As formulated, the parameters a, (3, have units of energy. The 
function U in (2.16) encompasses both the Schrodinger quantum stationary state and GL superconductivity models. 
In the former case, a(r) = V(r) — E, and (3 = 0, with E the energy eigenvalue and V the potential. For the GL 
model, a and (3 are independent of r, n* is the SC pair density at T = 0, and tp(r) is the complex order parameter. 

Since it is generally easier to solve the linear Schrodinger equation for the complex wave function than the coupled 
nonlinear Eqs. (2.8) and (2.9), we do not pursue the coupled approach here; however, there are some advantages to 
this approach for interacting electron systems 8,9 . The remainder of the paper concerns superconducting systems with 
a and (3 independent of the coordinate r. 

For mathematical expediency it is convenient to normalize the order parameter using a positive constant ip a . Since 
the GL energy and resulting GL variational equations are independent of the value oftp , its choice is arbitrary. Using 
the scaled temperature notation t = T/T c , where T c is the temperature at which a second order phase transition, 
with limit — > 0, would occur in a bulk sample, we employ the implicit temperature dependent normalization 
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The unitlcss temperature dependent function s(t) is defined by the fundamental length ratio 



s(t) 



where the lengths are defined as 
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(2.17) 



(2.18) 



(2.19) 



In conventional superconductivity, with (a < 0,(3 > 0), the length l a is referred to as the GL coherence length, 
usually denoted as £. For a bulk superconductor with A = the pair density for minimum energy \ip m {t = 0)| 2 = 
— a(0)/(3(0) = 1. Thus the maximum value of s(t) is s(0) = 1, i.e. l a (0) = 1/3(0). When A ^ 0, anomalous SC 
states exist with s(t) > 1, even when t > 0; thus we retain the notation l a to denote a fundamental length that can 
assume values outside the range of the conventional coherence length £. We also define a normalized, unitless, order 
parameter 



X(r) 



|V(r)| 



/n: 



s(t)f(r). 



(2.20) 



Because the scaling parameter n* is temperature independent, x is always physically meaningful; whereas / may 
become singular when the temperature dependent scaling parameter s(t) — ► 0. Thus the function /, although 
mathematically useful, is not in general a physical order parameter. 
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Since AG depends not only on the magnitudes, but also the signs of a and f3, it is convenient to write a = |a|sgn(a), 
and (5 = |/3|sgn(/J), where sgn(...) = ±1. Using Eqs. (2.16)- (2.19), the equations (2.8)-(2.10) in the volume V assume 
the form 

ll V 2 / - [sgn(a) + sgn(/3)/ 2 + l 2 a Q ■ Q]f, (2.21) 
V • (/ 2 Q) = 0, (2.22) 



s 



2 



V x [V x (Q - Q a )] + ( — ) f 2 Q = V x [V x {V(9 - a )}]. (2.23) 



with \l = (l/e*)y/m* / (jJ^n*) the London magnetic field penetration depth which is temperature independent. The 
inhomogeneous term in Eq. (2.23) is present when vortices, or 3D worm holes exist 9 . Equation (2.23) may be written 
in the form 

The right hand side of Eq. (2.24) is the physical super current density j, which we write in terms of a unitlcss current 
density parameter J by 

\ 2 

S 



j = -jj, where J = l a f 2 Q, j a = — — — — . (2.25) 

e*n l a \XlJ 

Although J is a convenient parameter for mathematical reasons, it is not appropriate for physical interpretation since 
it is singular in the limit s — > 0. Thus we also define the unitless current parameter 



For plotting, kj is a well behaved parameter that tracks the physical current density j. 
The surface boundary conditions from Eqs. (2.11-)(2.13) are 



l(a + W + (.-I,/<S M -£«^, ( , 27) 



where a = 1 + c — e, and 



h 2 

n x (H - H„) = — ^VaA, (2.28) 

Equation (2.22) is the continuity equation for j, and Eq. (2.29) is the boundary condition for the current normal to 
the surfaces. If A(f,9) contains a product of / and 9, then the boundary conditions (2.27) and (2.29) are coupled. 
From Eq. (2.14), the minimum free energy AG = G sc — G n for superconducting systems assumes the form 

AG mm = i^d 3 r[(A-A a ).j-sgn(/?K|/3| s 4 / 4 (r)] + ^ J ds ■ [nA + ^z< S 2 V/ 2 ]. (2.30) 

To proceed further one must formulate a model for the function A(/, 6), which determines the boundary conditions 
for ip and the minimum free energy. When currents are present a general form for A(/, 9) is not known. However, if 
the surface is that of a SNS junction, or Josephson junction, with tunneling supercurrent density j • n = j c sin A6, it 
follows from Eq. (2.29) that A must contain a term oc —j c cosA9, where A9 is the phase difference across the link. 

In the absence of transport currents, the superconducting sample in volume V is in thermal equilibrium with the 
surrounding material. One may assume that the sample and its surface are characterized by an order parameter %p, 
with A a function of \ip\ — ip f only. As discussed after Eq. (2.15), the surface energy term A contributes to AG mjn 
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only if A deviates from the form A oc \tp\ 2 . Thus, to obtain the leading surface contribution for J = 0, we introduce 
the surface energy function A as 

The first term is the form introduced by de Gennes 3 , where &2 is a characteristic length for the surface. The second 
term represents the superconducting pair interaction energy, with 64 a characteristic length. This term must be 
present in our approach for temperatures well below the critical temperature to be consistent with the retention of 
the / 4 term in the bulk free energy. The / 4 term may also be significant over most of the temperature range when 
the surface is strongly superconducting due to contact with a higher T c superconductor than the sample considered. 
Substituting the surface model (2.31) into Eq. (2.30), and using the boundary condition (2.27) to eliminate l/6 2 leads 
to the form 

AG min = - l -n* s \(i\s 2 |sgn(/3) S 2 J d 3 rf(r) + 1% J ds ■ f 

Although the volume term with a cancelled and the surface term with & 2 was eliminated, AG min remains an implicit 
function of a and 62 since / depends on them via Eqs. (2.21) and (2.27). If a 7^ 1 (c ^ e), it will be shown that 
AGmin is greater than it is for a = 1. Furthermore, only a — 1 permits a second order phase transition. For a = 1 
and ip/b± > 1 the parameter (3 can be negative and still have AG m „ < 0, i.e. the overall sample is in SC state. If 
i|/&4 <C 1, the surface effect is negligible and the parameter f3 must be positive for a SC state to exist. In the Gor'kov 
derivation 10 of the GL equations from the microscopic BCS theory 11 , (3 is always positive. However, surface effects 
were not considered in the Gor'kov derivation. 

As final observations for Section II, we examine the scaling roles of the fundamental lengths l a and lp, and the 
uniform solution of Eq. (2.21). It is evident from Eqs. (2.21), (2.23), and the definition of Q that l a scales the 
coordinate r, i.e. one may use r — ► r/l a . Multiplying Eq. (2.23) by Z 3 to scale the coordinate, the only parameter in 
Eq. (2.23) is the experimentally observable parameter k = X s /l a = Xl/1/3, where A s = \l/s is equal to the Ginsburg- 
Landau penetration depth Xql in conventional superconductivity. However, for the anomalous SC cases analyzed in 
section III, A s can lie outside the range of the conventional Xgl- The lengths l a and lp defined in Eq. (2.19) are 
positive real numbers. All signs consistent with the general GL Eq. (2.21) are in sgn(a) = ±1, and sgn(/3) = ±1, 
which dictate whether a solution / is an oscillatory or non-oscillatory function of the coordinate r. 

Consider the possibility of a uniform solution / u (r) = f of Eq. (2.21). Setting V 2 / = and using Eqs. (2.20) and 
(2.26) gives 

Q u = ±qq u , q u = -^ = I ^-[sgn(a) S 2 + sgn(/?) X 2 ] (2.33) 
From the definition of Q, Eq. (2.5), the phase difference over path r = to r = d is 

e{d) u -e{Q)u = -— dr-A±q u d. (2.34) 

H>o Jo 

Since the constant q u is real, the case with sgn(a) = sgn(/3) = 1 is excluded. Further restrictions are given in Eq. 
(3.16) below. Using Eq. (2.34), the uniform solution has the general form 

ip(r) u = exp j dr ' A ) [ c + etq " r + £-e~ lqur ] (2.35) 

As will be shown, the uniform solution plays an important role in the general solution space. This solution appears, 
not only at a second order phase transition where \tp\ — ► 0, but also well below the transition point where \ip\ is large. 

III. ONE DIMENSIONAL SYSTEMS 

Analytic Solution Classification 

In this section a complete set of analytic solutions of the first GL equation is derived for one dimensional systems, 
which in accordance with Eq. (2.22) have uniform current density. In general Eqs. (2.21) and (2.23) are coupled since 



\df f 
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(2.32) 
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Q remains multidimensional. However in SC microwires, 12 and certain thin films of thickness < Xql of arbitrary 
widths 13 , the current density is approximately uniform. Consider a one dimensional system with coordinate x. 
Approximating the magnitude of the true current density by its mean value over all x, and using J = l a f 2 Q, the first 
GL equation (2.21) is 



f a fJ I = sg n(a)f + sgn(P)f + J 2 f- 



subject to the boundary condition 



dx 2 



= -T-T s2f2 



6) 



(3.1) 



(3.2) 



where the subscript b denotes the surface boundary. At this point the rationale for introducing the normalized function 
/ and current parameter J is evident. Since one may use x/l a as a normalized coordinate, the only parameter in Eq. 
(3.1) is J. The unitlcss current parameter J, together with sgn(a) = ±1, and sgn(/3) = ±1 completely characterize 
the solutions of Eq. (3.1). 

When dealing with current carrying states, it is useful to introduce N{x) = f 2 (x). Equations (3.1) and (3.2) assume 
the forms 



Nli 



d 2 N 



2 J 2 + 2N 2 [sgn((3)N + sgn(a) 



la—\ 

dx ) 



and 



\{a + 1)N-^ + 2(o - 1)[J 2 + iV 2 (sgn(/?)A + sgn(a))] fl 2 a N^ 



1 1 

b 2 b A 



-^-^s 2 N\ 



(3.3) 



(3.4) 



which will be used below. 

Integrating Eq. (3.1), assuming J constant gives 



2 l a 



sgn(/?)/ 4 + 2sgn(a)/ 2 - 2J 2 /" 2 + C, 



(3.5) 



where C is a constant of integration. All solutions of Eq. (3.5) may be expressed in terms of Jacobian elliptic functions 
pq[u(x)|m], which are periodic in the argument u{x). In general when J ^ 0, determination of these functions requires 
the roots of a cubic equation in f 2 . However, the process is considerably simplified by dividing the solutions into two 
classes: 1. Functions with at least one positive finite extremum at x = xo, at which point C is evaluated. 2. Functions 
with f(x = x w ) = 0, at which point C is evaluated. In the first class only roots of a quadratic equation are required. 
For the second class only cases with J = are possible, and these may also be resolved by a quadratic. 



Solutions with Finite, Positive Extrema 

If / has at least one finite extremum point within each period, the constant C may be evaluated at an extremum 
point x 0} where df /dx = 0. Thus, we write Eq. (3.5) in the form 

2/2 2 = (/2 ~ ^){^(P)f 4 + M0)fo + 2sgn(a)]/ 2 + 2J 2 / " 2 }, (3.6) 

where f = /(x ) is an extremum to be determined from boundary conditions. A second extremum ff = f 2 (xi), 
given in terms of /o by setting the factor {...} = 0, may exist for both real and imaginary f\. 
Setting N(x) — f 2 (x), Eq. (3.6) assumes the form 

= 2 ( iV - Ar o){sgn(/3)7V 2 + [sgn(/3)A + 2sgn( a )]7V + 2J 2 A - 1 }, (3.7) 

Since the solutions are referenced to the extremum value N — N(xo), it follows from Eq. (3.3) that the sign of the 
parameter bo, given by 
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d N 2 

5 " = 'a-^-Uo - ^ 2 - Co], with Co = -iV 2 [sgn(/?)iVo + sgn(a)], 

determines whether No is a maximum or a minimum. 

For mathematical expediency, it is useful to make the following variable shift: 



Equation (3.7) yields 



which has roots 



with 



N(x) = No + sgn(b )v 2 (x). 



2 1 la ^) = s g n (Msgn(/3)(^ 2 - r+)(v 2 - r_), 



r± = 7 ± A, 



7 = -sgn(& ) 



-AT + sgn(a)sgn(/3) 



, A = V7 2 ~ &oSgn(/3) = A / — [d - J 2 sgn(/3)] , 



(3.8) 



(3.9) 



(3.10) 



(3.11) 



Co - -iV 2 [sgn(/3)iVo + sgn(a)], C\ 



N 



N 



+ sgn(a)sgn(/3) 



For J = 0, one must use the magnitude of C\ to obtain the correct result A > 0. The solutions of Eq. (3.10) are the 
various elliptic function combinations, listed in (3.13) and (3.14), whose specific form is dictated by the characteristics 
of the roots shown in (3.12). 

Roots: The array below shows the root ordering for all sign combinations that yield real N(x) > solutions. The 
dashed line separates the positive (5 cases from the negative (3 cases. 



Case 


sgn(a) 


sgn(/3) 


sgn(6 ) Root Order 


A^o Range 


J 2 Range 


A 


-1 


1 


-1 


r + > 0, r_ < 


< Ao < 1 


J 2 < Co < Ci 


B 


-1 


1 


1 


r + > r_ > 


0<N < 2/3 


C < J 2 < C l 


C 


-1 


1 


1 


r_ < r + < 


N > 2/3 


C < J 2 < Ci 


D± 


±1 


1 


1 


r_ = r* + 


< N < oo 


C < Ci < J 2 


E 


1 


1 


1 


r_ < r + < 


< N < oo 


J 2 < Ci 


F 


1 


-1 


1 


r+ > 0, r_ < 


< A^o < oo 


J 2 > Co 


G 


-1 


-1 


1 


r + > 0, r_ < 


< A^o < oo 


J 2 > C 


H 


1 


-1 


-1 


r + > r_ > 


A^o > 1 


J 2 <C 


I 


-1 


-1 


-1 


r + > r_ > 


< AT < oo 


J 2 <C 



(3.12) 



Solutions with J ^ 0: The function v{x) are Jacobian elliptic functions of the form pq[M(x)|m], shown in (3.13) 



v{x) 



r+\r-\ x-xo r+ + |r_| 

sd(w|m) — : — ' 



Cases 

A, F, G 

B, H, I 
C, E 



y|f+| sc(u|m) 
D± vV+| sc(w|m)dn(w|m) 



u(x) 



r + + |r_ 
JrZ sn(u|m) 



x-xo jr+_ 



x-x p I \r 

lev 



m 



x - xp \r+\ 
L f 2 



r+ + \r-\ 

r_ 

r+ 

\r+\ 
|r_| 

Re(r+) 



(3.13) 



1 - 
1 + 
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where \r±\ = —r± if r± < for real roots, and \r±\ = ^r + r_ = y/b~o for the complex root case D. The function 
N(x) = f 2 {x) is given by Eq. (3.9). Solutions for a conventional SC are cases A-C and D with sgn(a) = —1 and 
sgn(/3) = 1, which are the same as those given in Ref. 14. Solutions for the remaining 6 anomalous cases fall into the 
same solution classes, however 7, b , Co, C\ and A are different functions of No, due to different values of sgn(a) and 
sgn(/3). 

Solutions with J = 0: It follows from Eq. (3.6), with J = 0, that f(x) may be expressed directly in the form 
f{x) = f pq(u\m). 

m No Range sgn(6 ) 

N Q 




2 -No 

2 -Np 

No 
0.5N - 1 



0<N <1 -1 
1 < No < 2 1 
No > 2 1 



No-1 

yiVoTT 5^!L + I 0<iV <oo 1 
I\o + 1 



(3.14) 



u 1 x-x N I- N 

HE — V 1_ T T-O^iVo °~ Na - 1 1 

HI dn ^ 1<N <2 -1 
l a V 2 0.5iV 

H2 cn X -^V^~l ^ N > 2 -1 



' JVo-1 

I cn + i Ar 0<iV <oo -1 

l« iVn 



The zero current solutions in (3.14) may also be obtained from Eq. (3.13). The dashed line in (3.14) separates the 
sgn(/3) = 1 cases from the sgn(/3) = — 1 cases. The latter cases are possible only when superconductivity of the surface 
of the specimen is strongly enhanced through the boundary condition by introducing a significant quantity of pair 
interactions. For sgn(6 ) = +1, (—1), the function f = f(%o) 1S a minimum(maximum), respectively. 

We emphasize that / 2 = N is employed only for mathematical convenience. It is particularly useful in numerical 
computation when one needs the precise range of iVo which determines the parameter m in (3.14). However, the 
function \ = s .f 1S the normalized order parameter of physical interest, with \ip\ 2 = n*x 2 the actual superconducting 
pair density. Numerous treatises based on the GL model of superconductivity use the parameter / to illustrate the 
behavior of superconductors. In general, this is not a good parameter to use, except as an intermediate mathematical 
step. As seen in (3.14) the parameter / can become singular in cases C2, E, H2, and I, and is larger than unity in 
cases CI and HI. As shown below, / — > 00 when s — > 0; whereas x = s f remains finite and well behaved. Similar 
comments apply to the use of the current parameter kj in lieu of J. Thus it is x( or \ip\) and kj that should generally 
be used for final illustration of superconducting effects. 

It is interesting and useful to examine special solutions that appear at certain points in the (J,N ,s) parameter 
space. At the end of Section II we determined the complex order parameter for the uniform solution \u- Now we 
determine Xu from the boundary conditions. Setting dN/dx = d 2 N/dx 2 = for all values of x, Eq. (3.3) requires 
that J 2 — Co, i.e. b = 0. It follows, using Eq. (2.26), x — sV^V, and the boundary condition (3.2) with a = 1, that 
the (xu,s u ) coordinates of the uniform solution Xu are 

Xu = s u (0) = V-h/b 2 , s u (kj) = s„(0) A /-sgn(a)[sgn(/3) + (kj/xl) 2 }- (3.15) 

Since x/ s — V~N, the argument of the root in Eq. (3.15) is 1/A^„. For the uniform solution to exist for all x, the 
sgn(6 2 ) = -sgn(&4), and the order parameter Xu is independent of the current parameter kj; whereas s u depends on 
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kj. Because s is real and positive, Eq. (3.15) requires the [sgn(a), sgn(/3)] combinations and kj ranges shown in array 
(3.16). 



(3.16) 



Applying the uniform solution condition b = to the elliptic function solutions in (3.13) gives N(x) — Nq, with 
m = 0, at the boundaries of the solution case domains indicated in Eq. (3.16) by a vertical line. When kj ^ 0, the 
uniform solution point s u approaches infinity in both the large and small limits of \, corresponding, respectively, to 
weak and strong surface interaction limits. 

Since s(t) = y/\a(t)\/\(3(t)\ by definition, Eq. (3.15) gives the relationship &tt = t u between the volume parameter 
ratio a/(3 and the surfaces parameter ratio &4, b 2 as a function of the current parameter kj. For kj = 0, which requires 
that sgn(a) = — sgn(/3), it follows that 



Cases 


sgn(a) 


sgn(/3) 


kj Range 


A B, C 


-1 


1 


k 3 > 


F H 


1 


-1 


< kj < x 3 


Gjl 


-1 


-1 


kj > x 3 . 



a(t u ) 



(3.17) 



where t u is the normalized temperature at which s(t u ) — s u . For a bulk SC with kj = 0, the ratio a(0)//3(0) = — 1 
and x(0) = 1- In contrast, the ratio a(t u ) / (3{t u ) given by Eq. (3.17) is generally not equal to —1. 

The uniform solutions and other special solutions follow from the general solutions listed in Eqs. (3.13) and (3.14). 
For J = they are listed in Eq. (3.18) below. 



Case 

A, I 

E, F 
A|C1 
F|H1 

C1|C2, D_ 

H1|H2 

C2|E, H2|I 



m/fo 

' x — x 

—' 

x — x 



cos 



cosh 



sec 



sech 



1 
1 

X — Xq 

x - x 



nc 



b 

X — Xq 



x(o) 



m 



1 
1 




1 



N 



1 
1 

2 

2 

oo 



(3.18) 



The vertical line separating cases indicates that the special solution occurs at the boundary between the case solution 
domains. The solutions with Nq — ► are second order phase transition solutions which satisfy the linearized form of 
Eq. (3.1); those with Nq = 1 are the x independent uniform solutions; those with Nq = 2 satisfy the nonlinear form of 
Eq. (3.1); and those with — > oo have s — and represent the transition between Jacobian elliptic functions with 
and without inflection. 



Solutions for J = with f(x = x w ) — 

When J = 0, Eq. (3.5) assumes the form 

2 (l a ^\ = sgn(/3)/ 4 + 2sgn( a ).f 2 + C = sgn(/3)(/ 2 - r+)(.f 2 - r_), 



(3.19) 



where 
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r± = -sgn(a)sgn(/3) ± y/l - Csgn(/3). 



(3.20) 



Although the constant C is arbitrary, the only solutions of interest that were not included in solutions with finite 
extrema are those for which / is zero at some point x w , with df/dx\x w ^ 0. In such cases, the constant C is 



C = 2[l a 



d f(x w ) 
dx 



> 0. 



The root ordering and corresponding solutions are given in Eq. (3.22) and (3.23) below. 



Cases 
A' 
D' 
E' 

f\ r 



Case 


sgn(a) 


sgn(/3) Root Order 


C Range 


A' 


-1 


1 


r + > r_ > 


C < 1 


D' 


±1 


1 


r_ = r* + 


C>1 


E' 


1 


1 


r- < r+ < 


C < 1 


F', r 


±1 


-1 


r + > 0, r_ < 


< C < oo 




m 




u(x) 


m 


V 


frZ sn(ii m) 


x — x\ [r+ 

la V 2 


r_ 
r+ 



|r+| sc(u|m)dn(u|m) 
x/l^+l sc(u|m) 



^ - ^i \r+\ 
L V 2 



^e(r+) 
|r+| 



a; — xi / r 



r+ + |r 

where r± are given by Eq. (3.20), with \r± 



r +\ r -\ i / I \ x-x 1 r + + \r- 
sd(w|m) — : — • 



\r+\ 
|r_| 



r+ + |r_ 



(3.21) 



(3.22) 



(3.23) 



-r± if r± < for real roots, and |r±| 



C for case D'. In the 



application considered next, some of the functions in Eq. (3.23) are used to construct surface "pre- wetting" solutions 15 
which are zero at a point x\ near a surface. Since / > by definition, only that part of the elliptic function that is 
positive represents a physical solution. The slope of f 2 {x) at x = x\ is zero, thus establishing a smooth transition to 
the region with / = 0. 



IV. APPLICATION: PLANE SLAB SAMPLE 

The structure used in this section is a symmetric plane slab of thickness d, with the origin x = at its center. 
First we consider the zero current slab body solutions with iVo = N(0) — / 2 (0) ^ 0, then the corresponding current 
supporting solutions, followed by the currentless surface "pre- wetting" solutions with f(x w ) = 0. 



Currentless Slab Solutions 



For J = there are eight solutions listed in (3.14). Since N(0) may become singular, we replace it with the well 
behaved normalized order parameter x 2 (0) = s 2 N(0), and use lp = sl a . Applying the boundary condition (3.2) at 
x = 0.5c?, setting xq = 0, and defining u s = u(x = 0.56?), leads to the boundary condition 



ar! + -O_ + ^ x 2 (0.5rf)--^, 



(4.1) 



where 



fl+ = - 



dx I dx 2 dx 



0.5d 
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When a = 1 (c = e), it is seen in Eq. (4.1) that only the logarithmic derivative x~ 1 dx/dx is specified at the boundary. 
As mentioned, and as will be shown, a — 1 is the important case permitting second order phase transitions. 
With considerable manipulation of the Jacobian elliptic functions, one may derive the following forms for fl±. 



Case 




n± 




A 


On, 

Zu s 


l + l(l±l)[cn(2u 8 ) 


- dn(2u s )] 


sn(2u s ) 


CI 


Or,, 


1- l(l±l)[cn(2u a ) 


- dn(2u s )] 


sn(2u s ) 


C2,E 


On, 

ZU S 


dn(2u 8 ) + l(l±l)[l 

Z 


- cn(2u s )] 


sn(2u s ) 


F 


2u s 


cn(2« a ) + i(l±l)[l 


- dn(2u s )] 


sn(2u s ) 


HI 


2u s 


cn(2u s )-l(l±l)[l 


- dn(2u s )] 


sn(2u s ) 


H2,I 


2u s 


dn(2u s )- |(1 ± 1)[1 


- cn(2u s )] 


sn(2u s ) 



(4.2) 



The argument u s = u(x = 0.5d) and parameter m are easily transformed in (3.14) using x 2 (0) — s 2 N(0). The 
boundary conditions listed by Eqs. (4.1) contain the set {x(0),s(t),d/lp,d/b2,d/b4}. Specifying the three scaled 
lengths, the boundary equations give the normalized order parameter x(0) at the slab center as a function of s. In the 
boundary condition (4.1) the cases C and H in (3.13) split into subcases, defined in (3.14). The sgn(a) and sgn(/3) 
remain those of the parent case. The cases in (4.1) with sgn(6 ) = —1 have superconductivity at the surface reduced 
from that in the sample bulk; whereas cases with sgn(fe ) = 1 have enhanced surface superconductivity. 

An alternative form of the boundary condition Eq. (3.2) is obtained by integrating the first GL Eq. (3.1) across the 
slab. Noting that df(x)/dx = —df(—x)/dx and f(x) — f(—x), one obtains for the case a = 1 the relation 

_ + - s 2 f(Q.5d) = - j2 /(0 . 5d) J Q dx[ S gn(a)f(x) + sgn(/?)/ 3 + J 2 / -3 ]. (4.3) 

Since Eq. (4.3) is exact within the GL model, it is difficult to reconcile it with a similar, but distinctly different 
expression involving the density of states given by de Gennes 3 , who uses the same arguments to derive the GL 
equations. 

Up to this point, all relations follow directly or indirectly from the free energy (2.7) with the models for U, Eq. 
(2.16), and A, Eq. (2.31). No approximations, assumptions, or other conjectures have been made. To solve the 
boundary Eqs. (4.1) we assume that lp is approximately temperature independent. Comparing the surface energy 
density with the volume condensation energy density, one notes that l/b 2 (t) and \a(t)\ cx s 2 (t) are both coefficients 
which scale \ip\ 2 , and that I/64 and \[3\ both scale \ip\ A . Hence, by analogy, we conjecture that l/bz{t) may be a 
function of s 2 (t) and 64 is approximately independent of temperature. Since we have no model for 62 (0 we will 
neglect the temperature dependence and use b 2 (t) « fe 2 , where b 2 is a controllable external parameter. [Assuming 
that l/b^if) oc s 2 leads to unreasonable results.] 

The parameters b 2 , and 64 cannot be chosen arbitrarily. The free energy must be negative for a SC state to exist, 
which limits the range of 64, and the boundary condition at a second order phase transition limits the value and sign 
of b 2 . From Eq. (2.32), the normalized, minimum free energy for a symmetric slab is 

T min (s) = (2.r) 2 ^p = -\ (£j ( X 4 )sgn(/3) - ^x\0.5d) 2{a l) X 2 (0.5d)fi_, (4.4) 

where 

The normalizing term K* d is the total kinetic energy of a non-interacting gas of SC pairs in volume V = area x d, 
with wavelength d at t = 0. The entire slab is in a SC state providing that the parameter r min (s) < 0. It is clear 
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that if a = 1 (c = e) the inequality T m i n < can be satisfied when (3 > for all 64 > 0, or for 64 < in a restricted 
range. When (3 < then 64 must be positive, but not too large. For o^l the key factor determining whether a state 
is superconducting is the sgn(T mi „). States that exhibit a second order phase transitions (SOPT) arc defined by the 
zero limit entropy difference between S and N states, 

t ^ rain n . min ~ /, -\ 

hm — = — lim — = 0. (4.5) 

x^o dt dt x^o ds v ; 

It can be shown that all states with a > 1 are potentially superconducting states, and in restricted s(t) value domains 
states with a < 1 may be superconducting with a first order phase transition. However only states with a = 1 have 
second order phase transitions, as defined by Eq. (4.5). In this paper we focus on the SOPT states and states that 
are continuations of the SOPT states corresponding to t he case a — 1. 

When J = the uniform solutions are x« = s u = \/— 64/62, which may exist for cases A, CI, F, and HI. The 
normalized energy for a = 1 has the value 



1 A 4 X ' 



r 

oil. 
Z \ t>2 



d\ 2 , 2d 



(4.6) 



This exact analytical result for the uniform solution gives a simple check on the numerical results. 

The point where x( x ) — is interpreted as a second order phase transition point, and the corresponding s-value 
is obtained from Eq. (4.1) by letting x(0) approach zero. As can be seen from Eq. (3.14) cases A and I give rise to a 
solution with maximum order parameter, while solutions E and F with minimum order parameter x(0) in the center 
of the slab. Cases A and I, with sgn(o;) = — 1, lead to a characteristic equation for s = s in the limit x(0) — > 0. In 
this limit, from (3.18), x( x ) = x(0) cos(£x/d), and the characteristic equation (4.1) with a = 1 is 

C o tan(0.5(o) = ^ with Co = ^s Q . (4.7) 

02 1/3 

Similarly, for a = 1 the cases E and F, with sgn(a) = 1, and \{x) = x(0) cosh((x/d) gives 

C o tanh(0.5Co) =-■£-. (4.8) 

As will be shown below, if s is interpreted as a function of temperature, Eq. (4.7) gives rise to a reduction of t cs due 
to a pulling down of the order parameter at the surfaces of the slab, while Eq. (4.8) leads to an enhanced transition 
temperature t cs = T cs /T Cl as previously shown. 4 ' 15 

Fig. 1 shows the E solution for a — 1 and for a = 1 ± 0.001 near the transition point defined by x ~~ * 0. The 
solid curves are the square of the normalized order parameter x(0, s) = x(0) plotted as a function of the length ratio 
s(t) — l(3/l a (t), and the dashed curves are the energy parameter T m i n (s). Only the curve for a = 1 satisfies the 
entropy condition (4.5) for a second order phase transition. For a fixed set of parameters, the a = 1 states always 
have a lower value of T m i n (s) than states with a =/= 1 at the same value of s(t). 

As we go below the second order phase transition point the order parameter increases, but the corresponding s(t) 
value varies in different ways depending on the type of solution of the nonlinear equations, Eq. (4.1). As can be seen 
from Eq. (4.4), the minimum free energy T m i n can be positive or negative depending on sgn(/3) and the value of 64, 
and implicitly on 62 through x- In the conventional GL approach (3 > 0. 

Figs. 2 and 3 show the zero current solutions of the nonlinear Eqs. (4.1) for a = 1 (c = e) with the scaled order 
parameter at the slab center x(0) (solid lines), and the scaled minimum free energy Y m in (dashed lines), Eq. (4.4). In 
Figs. 2-10 the ratio d/\b2\ = 0.01. For a second order phase transition the sign of 62 is dictated by the linearized Eqs. 
(4.7) and (4.8). The value of d/lp = 0.5 in Figs. 2-10, i.e. the value of is fixed relative to the slab thickness d. 

Figure 2 shows our results for cases A, CI, C2, and E, all for (3 > 0, obtained from Eq. (4.1) for various values 
of b 2 and 64. The normalized energy function T min is shown by the broken lines and x(0), the normalized order 
parameter at the center of the slab, defined by Eq. (2.20), is depicted by solid lines. In this particular situation 
T min is always negative, but very small; thus negative values of 64 could drive the slab into the normal state with 
Trnm > 0. The regions over which the different solution types are valid are separated by a bar, and as can be seen 
from the intersection of both A and CI curves. Depending on the values of 62 and 64, the same type of solution can 
appear in different regions of the x(0) versus s diagram. Solutions A and CI meet when A^o = 1, at which point 
Xu = s u = \f— bi/b 2 = 1/V5 is a constant [See Eqs. (3.15)-(3.18)], independent of the coordinate x. As s varies 
through this point, the slope of x{ x — d/2) changes sign as it goes through zero. 
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The point where x(0) = is interpreted as the second order phase transition point and for a finite value of 62 it 
occurs at a positive, non-zero value of s. For solution A, the order parameter at the surface of the slab is "pulled 
down" with respect to the center of the slab while for solution E it is "pulled up" . We shall interpret below the A 
case as a decrease of the superconducting transition temperature T c due to the proximity effect caused, for example, 
by a normal metal deposited on the surface of a thin film, while the opposite applies to case E. The symbol for the 
reference temperature is T c , at which point x(0) = occurs when l/b 2 = 0. For case A, as s increases, x(0) = st/No 
increases, which is interpreted as a decrease in temperature, away from the phase transition point. 

For case E the value of s decreases as the order parameter \ increases. For s — > the value of x(0) remains finite 
when l/i>2 < 0. This happens because iVo = N(0) = ,/ 2 (0) — ► 00 when s — > 0. For s = solution type E and C2 
are satisfied by N — > 00 with Jacobian function parameter m = 0.5. The order parameter function x( x ) i s ! however, 
well behaved for s = as seen in (3.18) for the C2|E point. As x(0) increases, the function of type E switches to 
a C2 function and the s-value increases, as is the case for the A type solution, moving to lower temperatures, away 
from the phase transition point (lower A branch) . The same interpretation is now applied to the upper curve as we 
go from type C2 to CI to the upper A curve. Thus, if the phase transition point for the E type solution is interpreted 
as occurring at an enhanced transition temperature T cs > T c , then the point s = anchors the temperature T = T c . 

As one goes from the phase transition point along the E, C2, CI, and upper A branches, the parameter N — / 2 (0) 
goes from zero to infinity at s = 0, with x(0) finite, for the E branch. For the C2 section the function No is also 
infinite at s = 0. As s and x(0) increase, C2 changes into CI at N = 2. At this point the Jacobian elliptic function 
becomes a trigonometric function [See Eq. (3.18)] and then changes again into an elliptic function in the CI region. 
At the point where the CI joins the upper A branch, 60 oc d 2 x/dx 2 \ x =o changes sign, and the solution x{ x ) becomes 
independent of x at the intersection. In domains C2, and CI the net slope of x( x ) at the surface x = d/2 is positive, 
while in the upper A domain it is negative. 

The effect of changing the sign of the surface pair interaction parameter 64 is also shown. When 64 is negative 
it increases the net slope of x( x ) at the surface, enhancing superconductivity there by lowering the total energy of 
the slab. For solutions E, C2, and CI the net slope at the slab surface x = d/2 is always positive, thus aiding 
superconductivity in the slab, while type A gives rise to a negative slope at x = d/2 which has an inhibiting effect 
on the superconducting state in the slab. As mentioned, where solutions A and CI meet, the net slope of the order 
parameter is zero and x( x ) is constant, independent of x. 

Figure 3 depicts solutions corresponding to Fig. 2 with a = 1, but with (3 < for solutions F, HI, H2, and I. In order 
for these solutions to be superconducting the total minimized free energy difference r mi „, Eq. (4.4), must be negative. 
Equation (4.4) with a = 1 shows that this can be achieved only if 64 > with the surface term (<i/&4)x 4 (0.5cf) larger 
than the volume term 0.5(d/lp) 2 (x 4 } ■ As for (3 > cases, the surface parameter & 2 determines the point at which the 
second order phase transition occurs. Again, the results of Fig. 3 are exact solutions of our generalized GL equations 
and are completely independent of the temperature dependence of s(t). Although the lower curve for solution type I 
can be interpreted in the same way as the lower A-type curve in Fig. 2 was interpreted, the rest of Fig. 3 has some 
notable differences. The phase transition point x(0) = on branch F can be interpreted as an enhancement of the 
transition point above the reference temperature T c , defined by s{t = 1) = 0, because the order parameter at the 
surface is "pulled up" at x = d/2. Where F joins branch HI, the order parameter x = \J—bi/b2 = 1/5 is independent 
of x; thus the slope at the surface changes sign at this point and HI becomes a "pull down" case before the reference 
temperature at s — is reached. As the order parameter x(0) increases from zero, its value is finite at s — as it goes 
from F type to HI and H2 and then to a different upper branch of the I solution. Along the HI, H2 and I curves, the 
order parameter x( x = d/2) is always "pulled down" with respect to x(0). At the point x(0) = V2s, where HI and 
H2 join, xi x ) becomes a hyperbolic function, (3.18), and at s = 0, where H2 and I join, the value of N = f 2 (0) — ► 00. 
For the given parameter values, all of the solutions are superconducting, which depend on a positive 64 of the pair 
interaction parameter near the surface. 

Although the F solution is the only one with the order parameter "pulled up" at the surface, all solutions are 
superconducting. In particular, the HI and H2 solutions are "pull down" cases existing above our reference point T c 
defined by s(T = T c ) = 0. Thus enhanced superconductivity can exist for j3 < as long as the surface pair interaction 
parameter I/64 is large enough to overcome the unfavorable effect due to sgn(/3) = — 1 in the slab volume. This is the 
case for SNS or SS'S junctions where the N or S' regions are embedded by a strong superconductor S. 

The results depicted in Figs. 2 and 3, for the order parameter x( x )i defined by Eq. (2.20) as a function of the 
length ratio s(t), defined by Eq. (2.18), with both x{ x ) an d s(t) positive quantities by definition, are exact solutions 
of our generalized, nonlinear equations for a slab with kj = 0. They do not depend on, or show how, s(t) is related 
to the normalized temperature t = T/T c , which is not included explicitly in the GL equations. Thus we must ascribe 
a temperature relation to s(t) so that our model elucidates the arguments for an increase or decrease of the critical 
temperature. Since s > by definition, the temperature function related to the E, F, and H solution types must be 
different from that of the A, C, and I solution types. We introduce the relation 
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where the upper (lower) sign denotes the temperature model with t > 1 (t < 1), respectively, and apply it to the 
exact results shown in Figs. 2 and 3. When x(0) = 0, the fundamental length ratios s±(t) assume the values s + (t cs ) 
for t cs > 1, while S-(t cs ) for t cs < 1. 

At the reference temperature f = 1, the fundamental ratio s(l) — 0. For this model there are two temperature 
ranges < t < 1 and 1 < t < t cs , corresponding to the same value of s in the interval 1 > s > 0. For t cs < 1, an 
increasing s(t) value corresponds to an increasing x(0) value such as seen in the lower branches A and I in Figs. 2 
and 3. The explicit equation for s_(i cs ) is (4.7), and fixed parameters d/lp and d/62, from which t cs < 1 is obtained 
using Eq. (4.9). For t cs > 1, an increasing x(0) value corresponds to a decreasing s(t) value such as is depicted in 
Fig. 2 by the E branch and in Fig 3 by the F, HI, and H2 branches. In such cases the explicit equation for s+(t cs ) is 
Eq. (4.8), from which t cs > 1 is obtained using Eq. (4.9). 

With these t cs values, which are embedded in Eqs. (4.7) and (4.8), the results of Figs. 4 and 5 are calculated from 
the s(x(0)) values shown in Figs. 2 and 3 using Eq. (4.9). For the E, C2, CI, A curve in Fig. 4 and the F, HI, H2, I 
curve in Fig. 5, the x(0) curve is continuous and finite at t = 1, while A = / 2 (0) is infinite at this point. Again one 
sees that it is not the f(x) function which is the correct physical order parameter but the \{x) function. 

The A|C1 boundary point in Fig. 4 and the F|H1 boundary point in Fig. 5 are the uniform solution points, with 
corresponding temperature t u . Temperature t u may be found analytically using Eq. (3.15) with kj = 0, and Eq. (4.9). 
The result is 

where the upper (lower) sign denotes t u > 1 (t u < 1), respectively. Using Eq. (4.10), the temperature t u = \/2/3 in 
Fig. 4, and t u = ^13/12 in Fig. 5. 

The conservative and convincing results shown by %(0,i) lend credence to our results that boundary conditions 
imposed on small superconducting specimens may change drastically their behavior, including increasing the transition 
temperature. Changing the slope of the order parameter at the surface of thin films in a magnetic field parallel to the 
surface effects also the magnetic field at which superconductivity nucleates 16 . 

Although the temperature model, Eq. (4.9), is not unique, it appears to be in accord with experimental results 
described by curve A (lower curve) in Fig. 4 when the transition temperature of a thin superconducting film is depressed 
by bringing it in contact with normal metals (NSN). The upper curve in Fig. 4 is also a reasonable result for a weaker 
superconductor S' in contact with a stronger superconductor S. The model temperature behavior is different for the 
N-metal of a SNS junction: s oc \ff (See Section IV). 



Current Supporting Slab Solutions 

Next we consider the non-zero current solutions for a — 1, which evolve continuously from the solutions plotted 
in Figs. 2 and 3 as the current parameter kj oc j increases from zero. Using the solutions from Eq. (3.13) in the 
boundary condition (3.4), with a = 1, yields the x(0) versus s(t) plots shown in Figs. 6-10. A typical nonlinear 
boundary equation that we solved is that for cases A, F, and G. It is 

l_rf_&o 1 sn(u s )cn(u s ) d 2 _ _d_ 

2 1/3 VA N s dn 3 K) + 64 s_ b 2 > l4 - iij 

where 

N s = N(x — 0.5d) = No + -^^sd 2 (u s ), u s - \j^sVA, to = 1 (l + ^) . 

Figures 6 and 7 show the connectivity of all of the current carrying states A-E, with sgn(/3) = 1, for various values 
of the unitless current parameter kj. The bars denote the boundaries separating solution types. The parameter set 
(d/lp, d/&2, d/64), and the kj = curves are those depicted in Fig. 2. In Fig. 6 the temperature t > 1 is decreasing 
to t = 1 as s decreases to s = 0. It is clear that a second order phase transition, characterized by \ — 0, does not 
occur when kj 7^ 0. The upper curves in Fig. 6 show the physical order parameter decreasing for a fixed kj as s and 
the corresponding temperature increase; whereas the lower part of the curves behave unphysically. At a fixed value 
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of s(t), increasing the current parameter kj decreases the order parameter (starting from the kj = curve) until a 
critical, maximum current is reached, above which the order parameter does not exist. This critical point is given by 
an E(D+) solution for higher (lower) temperatures. In Fig. 7 the temperature t < 1 is decreasing from t = 1 as s 
increases from s = 0. Again, for a fixed s(t) there is a critical, maximum current parameter kj above which the order 
parameter does not exist. The critical point is given by a D_(A) solution for higher(lower) temperatures. In Fig. 7, 
the horizontal line \u — VT/5 is the value of the uniform solution, Eq. (3.15). The uniform solution is located at the 
A|C (A|B) boundary for higher (lower) temperatures. 

Figure 8 shows the A solution, which has a second order phase transition for kj = 0, for various values of kj. As s 
increases from s(t cs ) the temperature t < 1 decreases. Other observations are analogous to those for Fig. 7. 

Figures 9 and 10 show the connectivity of all of the current carrying states F, G, H, and I, with sgn(/3) = — 1, for 
various values of kj. The parameter set (d/lp, d/62, d/b±), and the kj = curves are those depicted in Fig. 3. In 
Figs. 9 and 10, the horizontal line is \u — 1/5, Eq. (3.15), which is located at the F|H and G|I boundaries. Other 
significant observations are analogous to those for Figs. 6 and 7. 



Simple Analytic Solutions 

Here we return to the first order differential equation (3.7) to develop a simple algebraic model that very accurately 
fits the kj = curves for the cases A, F, H, and I, which have two finite extrema. Denoting the second extrema by 
N\ = N(x\), where dN(x\)/dx = 0, the term in the bracket {...} in Eq. (3.7) is zero. Using x 2 = s 2 N and kj = s 3 J 
one obtains the exact relation between the extrema xo = x( x o) an d xi = x( x i)- It is 



[xl + xl + 2sgn( a )sgn(/3) S 2 ]x 2 ,X! + 2fc 2 sgn(/3) = 
Consider the kj = cases with sgn(a)sgn(/3) = —1. Equation (4.12) reduces to 



(4.12) 



Xo(s)+Xi(s) = 2s 2 , for XoXT^O. 



(4.13) 



It is evident that as s — > 0, the second extremum \\ — ~Xo- Since x 2 1S proportional to the physical electron pair 
density and is thus positive, Xi( s = 0) < is purely mathematical and does not represent a SC state. Assume x? 
may be expanded in a power series in s 2 . Taking the leading terms, we write 

Xl(s) -ao+ais 2 , (4.14) 



which requires that 



Xo(s) = -ao + (2 - ai)s 2 



(4.15) 



The expansion coefficients are obtained by evaluating Xo( s ) a t the SOPT point s = s Q , given by Eq. (4.7) or (4i 
and at the uniform solution xo = s u — \J — bi/b2, which gives 



a 



1- ( — 

So 



1 + 



1 - I — 

So 



(4.16) 



When s u < s a , the coefficient a < and a\ > 2. This is the F case, and Eq. (4.15) lies within plot accuracy for all 
s of the exact F-H1-H2 curve shown in Fig. 3. When s u > s a , the coefficient a > and a\ < 2. This is the A case, 
and Eq. (4.15) lies within plot accuracy for all s of the exact A-Cl SOPT curve shown in Fig. 2. Other curves that 
do not have a SOPT may also be obtained using Eq. (4.15); however the coefficients cannot be analytically evaluated. 

For kj ^ we have not found a simple expression for xi( s > %) snch that Eq. (4.12) produces the results shown in 
Figs. 8-10. It is evident from Fig. 10, for example, that an infinitesimal current completely changes the x(0, s) curve. 



Currentless Surface "Pre- Wetting" Solutions 



In this subsection we apply the solutions from Eq. (3.23), which exist only over a section of the slab. Substituting 
the solutions of Eq. (3.23) into the boundary condition (3.2), with a = 1 gives the equations to be solved for the 
integration constant C = 2{l a df / dx\x w ) 2 for a given value of the surface pre-wetting depth Iq = 0.5c? — x w . A typical 
equation is that for the E' case: 
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1 + ^/T~C dn(u s ) 



cn(u s )sn(u,s) 64 



+ -s 2 (l-VT~ C)sc 2 (u s ) 



d 



(4.17) 



where 



l 1 + 



m 



1 + VT~C' 



If the 64 surface interaction term is neglected, cases A' and E' are elliptic function representations of the extensive 
numerical solutions of Eq. (3.5), with J = 0, evaluated in Rcf. 15 by Indekeu and coworkers. 

For all cases, when the argument u(x) <C 1, the function f(x) oc u(x). In this limit, i.e. when the pre-wetting depth 

<^ la, the boundary condition (3.2) leads to ^ = —b 2 - 

Cases A', E', F', and F exhibit second order phase transitions (SOPT) in the limit C — > 0. Then the boundary 
condition for A' and F reduces to the form 



s (t) cot 



l s o (t) 



and for cases E' and F' it gives 



s (£)coth 



h 



l S o (t) 



d_ 



d_ 



for f(x w ) = 0, 



for f(x w ) = 0. 



(4.18) 



(4.19) 



It is interesting to compare the pre-wetting solution characteristic equations for lo = 0.5d with SOPT Eqs. (4.7) and 
(4.8) for the slab solution cases A, I and E, F. First note that in the A', F cases 62 < 0; whereas 62 > for the A, I 
cases. At the reference temperature, where s(t = 1) — » 0, Eq. (4.7) for cases A, I requires l/6 2 — > 0; whereas for A', 
F, Eq. (4.18) gives I&2I = lo at t = 1. Similar comments apply to cases E, F versus E', F'. 
The pre-wetting solution normalized, minimum free energy, corresponding to Eq. (4.4), is 




y [°' 5d dx X \x)sgn{P) - y X 4 (0.5d). 
h Jo.5d-l 04 



(4.20) 



As an application, consider a slab of normal, or SC, material between two superconducting reservoirs with a T c 
higher than that of the slab. In the absence of a transport current, the order parameter is characterized by an A type 
solution in the interior of the reservoir with extremum value /o < 1. At the surface there will be some bending such 
that f(d/2) < f . For the slab we assume an E' type pre-wetting solution with surface value 



f(d/2) = (l - VT~C) sc 2 (u s \m) 



(4.21) 



Using Eq. (4.21) with elliptic function identities to write the elliptic functions in the derivative boundary condition 
(4.17) in terms of C leads to the expression 



CM = f(d/2) 



2-/ 2 (d/2) 



(4.22) 



The value of f(d/2) is a function of s(t), but if one fixes the values of f(d/2) in Eq. (4.22), the resulting C(s) may be 
inserted in Eq. (4.21), to determine the pre-wetting depth Iq as a function of s(t). 

Fig. 11 shows the pre-wetting depth Iq plotted as a function of temperature, using t 2 = (1 — s 2 )/(l + s 2 ) for several 
values of f(d/2). Near T c of the SC reservoirs the surface value f(d/2) <C 1 and lo/d approaches \b2\/d — 0.2. As the 
temperature is decreased both f(d/2) and l have the tendency to increase until lo/d = 0.5, at which point we expect 
a transition to the E type. 



V. TRANSPORT CURRENTS IN MULTIPLE LAYERED SYSTEMS 



In the previous section we applied a boundary condition that was formulated from the surface energy characterized 
by the parameters hi and 64, which were assumed to represent not only the surface properties, but also the material 
beyond the surface. Now we formulate a set of boundary conditions that relate parameters across a boundary, 
bypassing the intermediate surface parameters. These boundary conditions are then applied to a SNS system with 
transport currents, and the results are compared with recent experiments. 
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Crossing Boundaries 



There are three basic quantities that are relevant at a boundary: the surface energy in G m i n , Eq. (2.15), the 
boundary condition (2.11) with a = 1, and the current density j, Eq. (2.24). The surface energy density is proportional 
to the function 



1 , ,,3 9 
9s = ~2^M dW\ 



A 



IVf 



1 

m* 



a 1 i n dA 



(5.1) 



and the boundary condition (2.11) may be written in the form 



-L[n.VhM 2 + 2A] =2g 8 . (5.2) 

Since g s is the same when viewed from either side of the surface at the same point on the surface, with m* changing 
value, in general, it follows that the right side of Eq. (5.2) must be continuous across the surface. Assuming that 
there is no generation or loss of particles at the surface, we conjecture that at the boundary between medium 1 and 
medium 2, the function A/m* satisfies 

A « - A < 2 > (5.3) 



m*(l) m*(2) 

Continuity of g s and A/m* gives the continuity of the SC density gradient condition 

M) Vxl{1)h - 5§ v * 2(2)l »- (5 - 4) 

The second boundary condition is the continuity of the physical current density j. Using the definition (2.26), 
continuity of j requires that 

tom-^MtMi). (5.5) 



m*(l)l p (l) Jy ' m*(2)l (2) 

Aside from known material and geometry parameters, the functions on each side of Eq. (5.4) depend only on s(t), kj 
and x( x o), where xo is the reference point in the medium. Thus, if j and t are fixed and the model s± is chosen for 
each medium, Eqs. (5.4) and (5.5) determine x(xo(l)) in one medium as a function of x( x o(2))- 

In order to determine the value of the pair {xo(l) 5 Xo(2)}, i.e. to set the SC level, a third constraint is required. 
If the system is open, the boundary condition (2.27) may applied at an outer surface; thus determining xo for the 
adjacent section as a function of the outer surface parameters. For a closed system, such as a ring constructed from 
two materials, the second constraint may be the flux quantization condition, which follows from Eq. (2.5), using Eqs. 
(2.25) and (2.26). It is 

-f'-(s?)=K" + >4 (5 - 6) 

where <f> is the flux enclosed by the contour I. Equations (5.4) - (5.6) determine the numerical value of x( x o) f° r each 
medium. 

In lieu of imposing an outer boundary derivative condition in terms of presently unknown characteristic lengths 
62,64, one may set the "level" by simulating the outer surface with the imposition of a discontinuity in the physical 
pair density \tp\ 2 at the inner boundary. That is, we set 

n:(l) X 2 (l)=r / n:(2) X 2 (2), (5.7) 

where the parameter 77 is a measure of the size of the jump in \ip\ 2 . Equations (5.4), (5.5), and (5.7) completely, and 
self-consistcntly determine the order parameter in each region of the system at a given value of kj and s(t) if the 
transition temperature of the system is known. 

Finally we note that the continuity of the current density j, also relates the gradients of the phases of tp(l) and 
■0(2) via 
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A Periodic SNS Layered System 



Here we analyze the experiments on SNS systems, AlAgAl and NbCuNb, as investigated by the Grenoble 
group 7 ' 17 ' 18 . Although both the (F, H) and the (E, D+) solution pairs are possible candidates for the N slab (Compare 
Figs. 6 and 9), we choose the (E, D+) solutions because they are pull up solutions throughout the entire temperature 
range, and because the positive value of (3 is expected to give a lower energy state. We believe that this choice is 
more likely to be the correct physical solution than that previously used 19 , which corresponds to the F solution. For 
simplicity, typical SC parameters and dimensions are used, similar to those used by Courtois et al. and Dubos et 
al. 7 ' 18 

We denote the superconducting layer (Al, Nb) properties by subscript "s" and those of the normal layer (Ag, Cu) 
by subscript "n" and use the Nq notation of Eq. (3.12) for the extremum of the f 2 function in the s and n layer. Note 
that for case A, the restriction is < No < 1, while for E and D + it is < No < oo. 

We assume that the current in the n layer kj n = — s^(t)J n is phase-coherent and is treated as a pair current 
as it is in the superconducting layers. We define in the n layer "coherence" length = l an = lj3 n /s n , with 

l\ n = (h/2it)(D /UbTc), where 3D = vp n £ n = ^ 2 /t, with D the diffusion constant, £ n the mean free path, vfu 
the Fermi velocity, and r the electron scatte ring time. The transition temperature of the SNS structure is T c , ks 
Boltzmann's constant and s n = \/i = y/T/T c . 

In the s layer the coherence length is £ s (i) = l as = lfSs{t)/ s s = £ s (0)(l + t 2 )/s s , with s s — V 1 — i 4 . The effective 
penetration depth A e // = ^IXl/V^- — t A , where Al is the London penetration depth and ft = 0.65(^o/Al) 1 / 3 > 1 for 
Al, a Pippard superconductor, for which £f(0) 3> £t)A| (£ = BCS coherence length) is satisfied. For Nb in the dirty 
limit il = (£ /6(0)) 1/2 > 1, with 

6(0) 



2nV2fi O H c (0)X eff (0y 

where (f> is the SC fluxoid quantum, H c (t) = H C (Q)(1 — t 2 ) is the thermodynamic critical field and no the vacuum 
permeability. The Drude resistance of the n layer R n = [m n /(n n e 2 T)](d n /A), with d n the thickness of the normal 
layer, A the cross-sectional area. 

The physical current density j is the same in the n and s regions, and thus the relation between j and the normalized 
current densities in the n and s regions is 

Ss \ 3 1 I fin T _ ( S s \ 3 1 I fin Al 



Jn UJ rip/ 8 [sj rUm + f 2 )V2H c (t)ns s h (5 ' 9) 

where r = (m s /m n )(n n /n s ). From the above, the physical critical current density j c , Eq. (2.26), leads to the energy 
expression 

ejcARn = ll{v) kjC = £ Wn{t) (5 - 10) 

where e = KD/d^ is the Thouless energy. 

It is our aim to calculate the largest (critical) current density kj c (t) = max[fcj] of a SNS structure as a function of 
temperature. The prefactor of kj C depends solely on properties of the n layer, but the unitless current parameter kj C 
tracks the physical critical current density j c and is implicitly connected to the SC and normal region parameters. 
The E and D + solutions applicable to the n region, and A solution relevant to the s region are, using Eq. (3.13) 

A N s (x)= N Q + ^sd 2 (u\m) 

E N n (x) - A + (l + ^iVo- A)sc>|m) t 5 ' 11 ) 
D + N n (x) = N + Vb~o~ sc 2 (u|m)dn 2 (u|m). 
For simplicity, we set m s /m n — n s /n n = 1 which reduce Eqs. (5.4) and (5.7), respectively, to 



dx 
and 



dN ni 

b — 



S s (t) 



S n (t) 
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2 

N s (b). (5.13) 

Although the parameter 77 is set by "external" boundary conditions, we assume continuity of the physical pair density 
across the boundary, i.e. \ip s \l = IVVil&j with 77 = 1. 

Since u and m are functions of Nq and J, it is possible to calculate for any pair (No, J) the values of N(x = d/2) — Nf, 
and dN/dx\(x = d/2) at the interface between the n and s regions. The triple (No, J, Nb) comprise a surface, and it 
is then possible to plot J vs. N for fixed Nb. 

The D + solutions for several values of Nf, for a NbCuNb specimen at T = 0.35 K are shown in Fig. 12. The D + 
solutions terminate on the right and the E solutions continue to J = [not shown] . For fixed No and Nb there exists 
a maximum J which is proportional to the critical current density j c (t) oc kj C (t), provided boundary conditions (5.5), 
(5.12), and (5.13) are satisfied. Only the curve with Nb = 16.73 satisfies all three boundary conditions at T = 0.35K ; 
thus the maximum value of this curve, which occurs in the D + solution region at T = 0.35i^T, corresponds to the 
unique j c . Similarly, Fig. 13 shows E solutions at T = 0.5 K for the above NbCuNb specimen. The E solutions 
terminate on the left and continue as D + solutions [not shown]. Only the Nb — 11.63 curve satisfies all boundary 
conditions, and therefore its maximum value, which lies in the E solution region, is the correct unique solution for the 
critical current at T = 0.5 K. The general behavior depicted in Figs. 12 and 13 is consistent with results expected 
from Fig. 6 for constant temperature. 

Critical current values kj c oc j c , Eq. (2.26), are plotted in Fig. 14 as a function of normalized temperature for the 
above NbCuNb and AlAgAl junctions with parameters listed in the captions for Figs. 12 and 14. The first three 
points on the left of each curve are D + solutions, and the remaining points arc E solutions. The plotted points are 
exact solutions with all three boundary conditions satisfied. The solid line is an aid to the eye only, connecting the 
computed solutions. The extrapolated values of kj c to T — K of the two junctions are 8.7 x 10~ 3 and 6.5 x 10~ 3 , and 
are comparable to the maximum kj value ~ 7 x 10~ 3 in Fig. 6, although they were obtained with different external 
parameters. Similarly, the maximum value of kj for the F type solution shown, in Fig. 9, is 8 x 10~ 3 . These results 
add credence to our earlier statements that it is kj oc j, and not J, that is a meaningful parameter. 

In Ref. 18 the dimensionless parameter eI c R n /e, Eq. (5.10), is plotted versus fc^T/e, with e a fitting parameter. 
These experiments indicate that this parameter extrapolated at T — K to the value 8.2. With the above parameters, 
the K limiting value of eI c R n /e for the NbCuNb junction is 47, and for the AlAgAl junction 257, indicating that 
this value depends strongly on the junction parameters of Eq. (5.10). 

Figure 15 shows the experimental points, taken from Fig. 4 of Ref. 7, of critical currents of a NbCuNb junction 
with parameters similar to those of Fig. 12. The solid line is calculated from the kj c (t) curve of Fig. 14 using Eq. 
(2.26) with e* = 2\e\, m*/n* s = 1.12 x 10~ 58 kgm 3 , and A = 5.9 x 10~ 14 m 2 . The latter curve is a convincing fit 
of our theory to the experimental points. Assuming that the extrapolation to T = K is meaningful, one obtains 
eI c R n /e = 21, using J c (0) = 2.5 mA, e = 24 /xeV, and R n = 0.2017. 

VI. CONCLUSIONS 

We have reformulated the GL theory by introducing the complete kinetic energy density, which requires a gradient 
term in the surface energy to support a SC state, thus replacing the standard GL energy density functional. The 
role of surface energy in determining the superconducting state of a sample was analyzed in detail. For a weakly 
superconducting surface, we have shown that the same phenomenological parameter 62 can be used to describe a 
reduction or an enhancement of the transition temperature in moderately small to small superconductors. Plating a 
superconducting specimen with a normal metal or ferromagnetic substance reduces the transition temperature and 
is described by 62 > 0. For 62 < an increase of T c is perhaps brought about by elastic strain 20 , observed on tin 
whiskers, by severely cold working the surface of InBi foils 6 , and by plating the specimen with a superconductor 
with a larger intrinsic T c than the specimen, and by other means 21 ~ 24 . In recent theoretical studies 25 the 3 Kelvin 
phase of Sr 2 Ru04 is modeled with Ru metal inclusion as interface states with locally enhanced (b 2 < 0) transition 
temperature. A negative bi value, or more precisely, a positive slope of the order parameter imposed at the surface 
of a normal slab, embedded between superconductors, induces superconductivity throughout the normal region if 
sufficiently thin. Nucleation fields for slabs for positive and negative slopes of the order parameter at the slab surfaces 
have been published previously 16 . They show that enhancements of H c ^ for 62 < is also possible. Other practical 
surface treatments and theoretical microscopic explanations relating to T c enhancements of superconductors are still 
to be discovered. 

The parameter (3, characterizing SC pair interactions in a sample, is always positive when the surface is weakly 
superconducting. If the surface is strongly superconducting, that is, when SC pair interactions play a significant role, 
superconductivity of the system is also characterized by a parameter b 4 . In this case, we have shown that the sample 



N n (b) = n 



s s (t) 



Sn(t) 
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may be in an anomalous superconducting state even when (3 is negative. The parameter 62 can be measured from 
SOPT experiments; whereas the parameter 64 can be determined from the uniform state, if it can be detected, or 
perhaps from the minimum surface free energy. 

For one dimensional systems with a uniform current density, all possible physical solutions of our generalized GL 
equations with sgn(a) = ±1 and sgn(/3) = ±1 were found and categorized. Although a commonly used parameter / 
is mathematically expedient, it can be misleading when used to depict the SC order parameter when surface effects 
are significant, i.e. / is not, in general, a physical order parameter. 

Introducing a transport current modifies the solutions of the theory with arbitrary parameters a and (3. Figures 
6 - 10, in which the order parameter x — ift/y/nt i s plotted as a function of the fundamental length ratio s = lp/l a , 
Eq. (2.18), relate to solutions with currents which are described by Eq. (3.14), while Eq. (3.23) and Fig. 11 are 
zero current pre- wetting solutions which evolve naturally in the present development. The current parameter kj, Eq. 
(2.26), tracks the physical current density. Since kj remains finite as s(t) — > 0, the critical current remains finite in a 
normal region, where s(t) = \ft, as t — > 0. In Figs. 12-15 the transport current results, applied to SNS junctions, are 
in excellent agreement with experiments over the entire wide temperature range of the measurements. The boundary 
condition for the gradient of the order parameter at the interface is derived from the surface energy functional. This 
boundary condition is consistent with that used in micronetworks 12 . The continuity of the order parameter at the 
interface, assumed without formal justification, completes the set of boundary conditions necessary to set the "level" 
of the order parameter. Although this condition is consistent with that used in quantum mechanics, the level of the 
order parameter is set, in general, by an external boundary condition on its gradient. 

It is often argued that the GL model equations apply only near T c . However, this argument is weak, since the GL 
phenomenological model does not contain the temperature explicitly. The exact solutions given here for the order 
parameter are functions of the fundamental length ratio s(t). It is true, if we write s = VI — t, as GL did, our 
results would be limited to temperatures near T c . However, there is no compelling requirement that restricts the 
temperature dependence of the present theory to this limited temperature range. The Gor'kov derivation 10 of the GL 
equations uses a small energy gap expansion of the microscopic BCS theory 11 , valid near T c , but that does not exclude 
the possibility that the phenomenological order parameter theory gives reasonable physical results well below 26 T c . 
However, it is prudent to relate s(t) to the experimental results of the temperature dependence l a (t) and lp{t). 



APPENDIX A: MINIMIZATION OF THE FUNCTIONAL G 



The energy functional Eq. (2.7) is 



G 



d v 



2^[VM • V|V| + |V>| 2 Q • Q] + ^(V x A') 2 



^ ^ds-[nA + (c-e)|V>|V|V>|]. 



where A' = A — A a . 

The total variation of G is 



2m* 



SG 



in 



dU 



v h 2 9|VI 
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HQ • Q]*M + V|VI • V5|V| + \M 2 Q ■ - 



e s — IVf Q • SA + — 2 (V x A') • (V x 6 A)} + 



L 



ds- [nM + (c-e)J(|V|VM)]. 



Consider the vector identities 



(Al) 



(A2) 



V/ • Vg - -/V • Vg + V • (fWg) 
B • (V x A) = A • (V x A) + V • [A x B] 



(A3) 
(A4) 



Using the vector identities for the volume terms involving VS\ip\,V89, V x SA, and applying the divergence theorem, 
one obtains from Eq. (A2) the variational form 
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9m* f m* 811 

-^SG = 2 J v d 3 r{\^Tg^ + \m ■ Q - V 2 |V|]^|V| - V • (H\ 2 Q)S9 + 
2ir m* 

I— m 2 Q+ —2 V x (V x A')] • 5 A} + (A5) 

J ds ■ {[(a + i)vivi + (a - + + m 2 Q + + 

2m* 

n[ j (V x A') x n + V A A] • 5A}, 

where a = 1 + c — e. The variations <J|^>|, 50, 5A may be taken independently in V and on S. Thus the coefficients of 
the variations are zero and one obtains Eqs. (2.8)-(2.13). 

To obtain the form of the minimum free energy given by Eq. (2.14), use (A3) with g replaced by / to eliminate 
V/ • V/, and use Eq. (2.8) to eliminate V 2 |# 
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FIG. 1. The order parameter x(0) of the E solution for is plotted as a function of the fundamental length ratio s(t) = lp/l a (t) 
for a = 1 and for a = 1 ±0.001, where a is a measure of the amount of Vx 2 in the surface energy. The dashed curves are the 
energy parameter Y min , Eq. (4.4). Only the a — 1 curve satisfies the entropy condition (4.5) for a second order phase transition 
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FIG. 2. The order parameter x(0) is plotted as a function of s for positive f3 cases A, CI, C2, and E for three surface 
parameter sets (62,64)- Also shown by dashed curves are the normalized, minimum free energy Y min . A superconducting state 
requires F min < 0. 

FIG. 3. The order parameter x(0) is plotted as a function of s for the negative /3 cases F, HI, H2, and I for two surface 
parameter sets (62,64)- The dashed lines are the minimized energy P m i n , Eq. (4.4). The superconducting state requires 

FIG. 4. The order parameter x(0) for the cases shown in Fig. 2 are plotted as a function of the normalized temperature 
t = T/T c using Eq. (4.9). 



FIG. 5. The order parameter x(0) for the cases shown in Fig. 3 are plotted as a function of the normalized temperature 
t = T/T c using Eq. (4.9). 



FIG. 6. The order parameter x(0) is plotted as a function of s for positive (3 cases E and D+ for different values of the 
current parameter kj. The kj — curves is that in Fig. 2, curve E. 

FIG. 7. The order parameter x(0) is plotted as a function of s for positive (3 cases A, B, C, and D_ for different values of 
the current parameter kj. The kj — curve is that in Fig. 2, upper branch C2, CI, A. 

FIG. 8. The order parameter x(0) is plotted as a function of s for positive (3 case A (lower branch) for different values of the 
current parameter kj. The kj — curve is the lower branch in Fig. 2. 



FIG. 9. The order parameter x(0) is plotted as a function of s for positive (3 cases F and H for different values of the current 
parameter kj. The kj = curve is the F, HI, H2 curve in Fig. 3. 

FIG. 10. The order parameter x(0) is plotted as a function of s for positive (3 cases G and I for different values of the current 
parameter kj. The kj = curve is the upper I branch in Fig. 3. 



FIG. 11. Shown is the pre- wetting depth lo of the E' solution, Eq. (3.23), plotted as a function of temperature for fixed 
values of f(d/2) at the NS interface of a SC reservoir. 

FIG. 12. Normalized current J vs. the minimum of the squared order parameter f§ — No in the n region at constant 
temperature for various N b — f£ values at the SN interface for the D + solution of NbCu at T — 0.35 K. The continuation of 
these curves on the right are E solutions [not shown]. Parameters: d n — 800 nm, d a = 200 nm, lp n = 61.3 nm, l/3 a (0) = 22.1nm, 
l„ = 48 nm, \ L = 39 nm, £ = 38 nm, T c = 8.1 K, v Fn = 1.57 x 10 6 m/s. 

FIG. 13. Similar to Fig. 12, but for the E solution of NbCu at T = 0.5 K. The continuation of these curves on the left are 
D + solutions [not shown]. The parameters are the same as used in Fig. 12. 

FIG. 14. Normalized critical current parameter kj C oc j c , Eq. (2.26), plotted as a function of normalized temperature for 
a NbCuNb junction with the parameters of Fig. 12, and the same for an AlAgAl junction with parameters: d n = 1.48 /im, 
d 3 = 2d n , lf} n = 113 nm, l fjs (0) = 487 nm, l n = 33 nm, X L = 14.8 nm, £ = 1-6 /an, T c = 1.45 K, v Fn = 1.38 x 10 6 m/s. 

FIG. 15. Critical current I c plotted as a function of temperature T of a NbCuNb junction. The experimental points are from 
Ref. 7, and the solid line is calculated from k jc of Fig. 14. 
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